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Abstract. In this paper we use spontaneous flux production in annular 
superconductors to shed light on the Kibble-Zurek scenario. In particular, we 
examine the effects of finite size, and of external fields, neither of which is directly 

amenable to the KZ analysis. We sec, from extensive simulations, that a realistic 2D 
superconducting annulus in 3D space, with dynamical electromagnetic fields behaves 
very like, and often indistinguishably from, a ID ring with classical EM fields. 
Further, the properties of that ID ring are well represented by analytic Gaussian 
approximations. Experimental results for anmili in the presence of external fields 
confirm these numerical findings, which have implications for the relevance of causal 
horizons to determine defect production. 
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1. Introduction 

Causality imposes constraints on systems that are strongly out of equilibrium by 
restricting the rate of change of correlation lengths to the relevant causal speed (such as 
the speed of sound). This is particularly the case in the vicinity of symmetry- breaking 
continuous transitions where, adiabatically, correlation lengths become arbitrarily large. 
Since such transitions take place in finite time, correlation lengths must remain finite. 
The faster the transition rate and the shorter the time available for system ordering, 
the more constrained the correlation lengths will be. 
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Figure 1. 3D sketch (not to scale) of a superconducting ring of width w with n trapped 
magnetic flux quanta, in the absence of any external flux. The dashed line indicates 
that a supercurrent, Idr, circulates around the ring. The magnetic field (solid) lines 
wrap around the ring section. 



The suggestion that causahty would constrain correlation lengths was originally 
made by Kibble [H [2] in the context of the very early universe. The relevance for the 
early universe is that if the system is frustrated this is relieved by the spontaneous 
creation of topological defects, whose separation reflects the correlation lengths at the 
time of their appearance, and which should be visible today. Zurek[3lll] independently 
reached the same conclusions for condensed matter systems in what is known as the 
Kibble-Zurek (KZ) scenario. 

This paper focuses on the behaviour of superconducting annuli. We follow Zurek in 
considering spontaneous flux production arising from rapid temperature quenches into 
the superconducting phase as a result of constrained correlation lengths. An idealisation 
of our experimental setup - which we shall discuss later - is shown in Figjl] Further 
details are given in the appendix. 

To model the system we assume a simple complex scalar field (pix) that can act as 
a proxy for a Cooper pair, decomposed as 

0(x) = |0(x)|e^^(-). (1) 

For the sake of argument we restrict ourselves to the inner radius of the annulus, where x 
measures the distance along the ring, circumference C. On quenching from the normal to 
the superconducting phase, we can define a winding number density n{x) = dxO{x)/2TT 
along this circumference. In the absence of any external fields the total winding number 
n around the loop is zero on average, but will have non-zero variance 

rC ,c 

(n^) = / dx I dy {n{x)n{y)), (2) 



Jo JO 

and it is this that we measure, indirectly, through the matching fiux generated in the 
interior of the annulus, which is directly observable. This fiux is quantised (as fiuxoids) 
in units of $o = hc/2e. 

For large annuli, for which C ^ f , we can replace ^ by 

(ji^) = j dx {n (x) n (0)) ^ — dx {n (x) n (0)) , (3) 
2 Jo 2 Jo 

if we assume that the two-point correlation function {n{x)n{y)) is of short range 

compared to C. Suppose the correlator is described through the single length scale 
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^. Simple dimensional analysis then gives 

{n') = a J (4) 

for ^ of ([s]). This is the perimeter law assumed in the Kibble- Zurek picture, interpreted 
as a random walk in phase along the circumference in steps of length ^. The same 
analysis would apply to spontaneous superflow in an annular tube of superfluid '^He 
[4J but whereas it is difficult to perform such a measurement there, in our case we can 
measure the flux directly to determine the characteristic correlation length. 

To estimate the length scale ^ we repeat the Kibble-Zurek argument of [H H], 
assuming rapid cooling through a continuous transition, transition temperature Tc. 
Suppose ^(t) is the adiabatic correlation length for n{x) at time t, diverging at the 
transition at t = as the temperature T{t) changes (T(0) = Tc). If the system is 
initially homogeneous and isotropic the KZ scenario proposes that on approaching the 
transition the increasing correlation length, originally changing adiabatically, will freeze 
(the 'impulse' regime) at time t< < 0, when ^(t<) ~ c(t<), where c(t) is the relevant 
causal speed at time t. After passing through the transition the system will unfreeze 
into a new adiabatic regime at a time t> > when |C(^>)| ~ c(t>). The correlation 
length ^ that sets the scale for phase change along the ring is 

e^e(t>) = eo(^y, (5) 

where and tq are system-dependent and tq is the quench time (the inverse quench rate 
through the transition). The exponent a can be derived from critical exponents valid 
for the adiabatic transition and the out-of-equilibrium behaviour is thereby categorised 
by the equilibrium universality class of the system. For Type-II superconductors 
t> = y/TQTQ and o" = 1/4 in the mean-field approximation. 

There is a somewhat different viewpoint, predicated on the approximate 
Gaussianity of fast quenches, explored several years ago by one of us (RR) [3 E] . More 
recently, preliminary simulations which further support underlying Gaussian behaviour 
were performed by us [H |8] and this paper continues this analysis in two ways. Firstly, 
we extend both the numerical simulations and the analytical approximations to see 
more closely the validity of simple Gaussian assumptions. Secondly we present new 
experimental data on the presence of fluxoid production in the presence of an external 
field, which go beyond the reach of the KZ scenario, but for which Gaussian assumptions 
give predictions. 

Our starting point is the observation that we would get the same scaling behaviour 
(|5| if we evaluated ^ at t< and, indeed, it has been proposed that t< is as relevant 
a time as t> [21 E], but we know for other reasons that it is what happens after the 
transition is implemented that determines the relevant distance scales [9]. This suggests 
defect formation as due to the growth of unstable long wavelength modes as the field is 
rolling off the hill which now denotes the unstable symmetry preserving vacuum, rather 
than through direct causal bounds. If, nonetheless, we take the KZ scenario as a guide 
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for the time t> at which defects form, then t> is much smaller than the time it takes 
for the order parameter field(s) to experience the degenerate vacua. If they remain 
so close to the unstable 'vacuum' that the non-linearities of the backreaction can be 
approximately ignored then the effectively linearised field equations give rise to Gaussian 
field correlation functions. Generally, these lead to the identical scaling behaviour (1) for 
fast quenches for large systems, both non-relativistic and relativistic [5l|6]. In practice, 
whatever the difference in language, our viewpoint is not very dissimilar from the KZ 
picture, in that both our simulations and the many simulations (e.g. [lOl IHl [T2] ) 
demonstrating the validity of the scaling behaviour ([S]) are based on time- dependent 
Ginzburg-Landau (TDGL) theory. It is these equations that we linearise to derive the 
Gaussian results. 

Having said that, there are two issues which force us to step outside the KZ picture 
and it is here that the greater generality of our approach is crucial. The first is that, 
whereas the early universe can be taken as homogeneous to a very high degree, the 
same cannot be said about laboratory systems. One step to limit inhomogeneity is to 
make the systems small. However, this looks to be swapping one problem for another 
in that we have now incurred finite size effects. Specifically, the scaling behaviour (|5]) 
assumed that C ^ ^. The Gaussian approximation permits an analytic replacement 
to Q and ^ for homogeneous small systems for which C < ^ which is inaccessible to 
the KZ picture. Thus, what seems to be a problem is, in fact, an opportunity. We can 
use small annuli to cast the KZ scenario in a different light and illuminate the Gaussian 
approximation. 

Secondly, the Gaussian approximation has a further reach, to explicit symmetry 
breaking driven by the application of external fields, for which the KZ scenario again 
offers little help. This is an interesting problem in its own right and experimentally 
is a serious problem because of the presence of stray fields in the laboratory. Scaling 
behaviour (|5| can only be demonstrated once external fields have been countered upon 
which, again, the Gaussian approximation gives insight. 

This paper is organised as follows. For narrow rings we might hope that the 
system behaves one-dimensionally. In the next section we reexamine the Gaussian 
approximation (or rather, two closely related Gaussian approximations) both for large 
and small idealised ID systems. We shall show later that numerical simulations provide 
strong support for Gaussian behaviour, both in the power spectrum of field ordering 
and in the damping of fluxoid generation in small systems, an exercise begun in [TJ E] 
and developed considerably here. In particular, we explore the way the power of field 
fluctuations is driven into long wavelengths through the growth of unstable modes, in 
good accord with our Gaussian expectations. 

In the subsequent sections we briefly summarise earlier work of ours [H] in which we 
simulated a 2D superconducting annulus embedded in 3D space, with a fully dynamical 
electromagnetic field in three dimensions. The reader is referred to that paper for further 
details. Surprisingly, we found that our results are, in key ways, indistinguishable from 
those of the ID simulation and thereby, of the Gaussian approximation. 
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We conclude with a discussion of new experimental results for fluxoid production 
in an external field and show that it too, is explicable in terms of Gaussian fluctuations. 
As we said earlier, the experimental details are given in the appendix. 



2. Gaussian behaviour for ID systems 

Unsurprisingly, ID systems are the easiest for which to make analytic predictions and to 
implement numerical simulations. There are two strongly related variants of Gaussian 
approximation that we shall consider; each has its strengths. 



2.1. Long superconductors: Gaussian probabilities 

Our first Gaussian approximation, first proposed in [131 [E], has the advantage of being 
much simpler analytically, predicated on the perimeter rule Q. We begin by considering 
a large loop, divided up into domains, in each of which 6' is a constant. We assume that 
there is no correlation between the values of 6 in adjacent domains but, in calculating 
the phase change A6 for the loop, the geodesic rule is adopted. This means that the 
shortest path in phase will be taken when jumping from one domain to the other. 
Let On {AG) be the probability that the change in phase 6 is after the N links that 
make the loop. The assumed lack of correlation means that the probability of ending 
with a phase shift of 27rm (i.e., net fluxoid number m) is: 

/7r+2r7i7r 
dQGNiO). (6) 
-7r+2m7r 

It follows that Pm{N) = p^rn{N). 

To bring this into correspondence with the KZ scenario, we should identify the 
domain size as comparable to C/^, i.e., C = aNC,, where a = 0(1). With this in mind, 
we assume that the total phase change A9 around the loop can be expressed as the 
sum of a random term and a geodesic-rule correction 6Q, necessary to obtain an 
integer winding number [15]. It is convenient to relax to be a continuous variable. 
We further assume that G has a normal distribution with average G = and variance 
cr^ proportional to A^ i.e.. 

The trapping probabilities pm = P-m are easily found: In particular, since in our 
experiments we rarely see more than one fluxoid, we are primarily interested in 

Po{N)=eTf\ , 1 (8) 

r Stt tt "I 

Pi = pi{N) +p-i{N) = erf , - erf , . (9) 



and 
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From the large-N behaviour of the pm it follows that, for large rings, cr^(A^) = 
47r^(n^). Also, we can identify Eq. ([t]) from the central limit distribution by taking 
cr^(iV) = for large A^. Further details may be found in Ref. |15j . 

2.2. A Gaussian variant: Gaussian distributions 

We stay with large annuli but, rather than assume the Gaussian distribution for B of 
([7]) we assume that the winding number density field n(x) is a Gaussian field. It then 
follows that the all important probabilities of finding no fluxoids or one fluxoid are 

fo = ^ r dz eM-^'{n')/2) (10) 

and 

/i = f, + f_, = -[ dz exp{-z^n'')/2) cosz (11) 

respectively. We use the notation fm for the probability of winding number m in this 
approximation to distinguish it from the Pm of the previous section for Gaussian winding 
number densities. Neither ([s]) or (|9]) take periodicity into account. 

In Fig. |2] we show how the predictions for p^ and /„ are almost indistinguishable 
once we get away from very small probabilities where periodicity is important. Further, 
we see that, with that caveat, they agree with the simulations for the fully interacting 
theory, discussed below. In particular, the Gaussian upper bound on finding a single 
fluxoid at about 0.49 that follows from ([s]) or ([9]) is saturated in the simulation (within 
errors). This supports the Gaussian approximation in its simplest guise. 

The Gaussian approximations so far do not give the KZ scaling behaviour of ([s]). To 
see this we make the further Gaussian approximation that, on decomposing the complex 
order parameter field as = (0i +202)/v^, 0i and 02 are independent Gaussian fields. 
The correlation function for the winding number density is now determined by the 
correlation function G{x) for the field components, defined by 

{(Pa{x)(Pt{y)) = SabG{\x-y\). (12) 

We now use the fact that, for damped systems, the Gaussian approximation 
corresponds to a linearisation of the time-dependent Ginzburg-Landau (TDGL) 
equations. Periodicity is unimportant for large rings and we find that, for the linear 
TDGL equations (TU], field correlation functions G{r) at the defect formation time for 
large systems take the form (r = |x|) 



G{r,t) = j 4ke'^-''P{k,t) (13) 
in which the power spectrum P{k, t) has a representation in terms of the Schwinger 



proper-time r (in the dimensionless units of (26) below) as (approximately) [5| l6l [T6] 



P{k,t)oc dr e-^''\- (14) 
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Figure 2. Wc plot (po, /o) and then (pi, /i) against (n^) for a ID superconducting 
ring to slrow tlie similarity of the Gaussian approximations. In particular, both pi 
and /i have maximum values of approximately 0.49. Both are contrasted with the 



results of the ID simulation derived from (26 1. Because neither approximation takes 



periodicity into account the likelihood of seeing flux at low levels is overestimated; see 
Fig. [5] Otherwise agreement is good within errors. 



where e(t) = T{t)/Tc is the reduced temperature. We assume the dependence e{t) = 
—t/Tq for a quench hnear in time in the vicinity of the phase transition where fluxoid 
formation takes place. Adopting this e(t) gives 



P{k,t) oc / dre 



V4T-Q 



(15) 



P{k) (X k ^ for large k"^. In the dimensionless units above, the time for the formation 



of fluxoids is ty = 0{y/fQ) 



There is a note of caution. ID systems have some peculiarities (such as the fact 
that formally there are no transitions in ID). In D spatial dimensions the integration 



variable ^k of (13) is replaced by ^^k. When D > 2 the integral is dominated by the 
ultraviolet for slow quenches, which needs to be cut off at distance scale ^o- 

To convert the Gaussian behaviour of the fields into the Gaussian behaviour of the 
winding number density is a messy process. A simple approximation, sufficient for our 
purposes, gives [7j 

{n (x) n (0)) ^ ^ [/' {xf - f" (x) f (x)] (16) 



(2vr) 



where / (x) 
For D 



--G{x)/G{0). 

1 or provided quenches are fast enough that we can ignore the ultraviolet 



effects then, evaluated at the KZ time t> 
/(r) ^ e-"'/2«' 



^1/2 



/(r) (r = |x|) takes the Gaussian form 

(17) 
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for small r, where ~ ^, the KZ separation length |5l E]. This is appropriate for a 
damped system at short distance. 

This creation of winding number can be thought of as due to the unstable long 
wavelength modes of the field growing as fast as possible, with no backreaction in the 
short time necessary. It follows that 

{n^)^a^r-^y\ (18) 

where a = 1/4. 

All the above ignores periodicity. The advantage of choosing Gaussian fields is that 
they permit periodic boundary conditions in a way that the Gaussian distributions as 
posed above do not. 



2.3. Small annuli: Gaussian approximation 

We would expect that small systems, whether superconductors or not, would show 
proportionately less defect production per unit length than larger systems because of 
end effects or periodicity. Let us consider small annuli for which C/C, ^ 1. 

The periodicity of f{x) in x (mod C) is now important. The effect of periodicity is 
to discretise in (15). To implement this in the approximation (17) we need to replace 
f{x) by its periodic generalisation, the Jacobi ^ function |7j 

_ M^rxlC]2n^aC^ .... 



^ 1 -4sm\nx/C) e'^- , (20) 

whence, from ([2]), 

{n') = 0(e-^-'«'/^'). (21) 

Thus, rather than the power falloff for large loops we have exponential damping 

In(n^) ^ - Att^C/C^ + const. (22) 

Since it is small, (n^) ^ fi, the probability of finding single winding number along the 
annulus. 

In principle C,r oc Tq, but some caution is necessary in that the approximation of 
(16) is oversimple, which makes the single term (22) only approximate. In we showed 
that exponential damping does take place. However, we can go further in that, if we 
combine (18) for cr = 1/4 with (22), it follows that 

{n') = F{tq/C') (23) 

for some universal function F interpolating between power behaviour and exponential 
damping that covers all sized annuli and quench rates. 
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2.4- Numerical simulations 

We simulate our ID superconductor with the U(l) scalar field theory of the complex 
order parameter 0, given previously, on a ring of circumference C with periodic boundary 
conditions. There is no way to include a dynamical electromagnetic field in a ID 
simulation. To mimic a temperature quench through its critical point at time t = 
we take an explicitly time-dependent potential 

Vm')=e{tM\' + h\<f>\\ (24) 



Rather than just take e linear in t, as we did in (15), we now adopt the more realistic 
behaviour 

[ 1, ^ < -TQ 

e{t) = I -t/TQ, _T-Q < t < tq , (25) 

[ -1, ^ > tq 

to model a slow quench, in which the quench has a beginning and an end. We start at 
t = —2tq and continue until t = Atq, by which time the defects have frozen out. This 
system is modelled by a damped second-order Langevin equation, 

dV 

- dl(t>a + VdfPa + TT" = Ca, (26) 
OCpa 

with zero mean Gaussian noise for which 

(Ca(x', t')Ux, t)) = 2r]T6{x' - x)6{t' - t)6ab. (27) 

The interaction strength b is taken to be unity and we take rj = 1 and T = 0.01, yielding 
a heavily damped system.. Our linearised Gaussian approximations are independent of 
b in the first instance. A mean-field generalisation of the Gaussian approximation in 
which the self-interactions are also linearised does give a weak 6-dependence (e.g. see 
[3]). We have performed our simulations for larger values of b without any noticeably 
different results. In evolving the equations a stochastic leapfrog method is used [IT]. 
Further details may be found in Ref. [7]. 

Our first application of this simple simulation in the current work is to test the 
validity of our Gaussian approximations discussed earlier. In Fig. |2]we show the observed 
frequencies for trapping a given number of fiuxoids as a function of (n^). This is 
compared with the Gaussian predictions {po, /o) and (pi, /i). Agreement is surprisingly 
good - at least for sufficiently large {n"^)- 

In Fig. [5] we compare the simulation results, scaled according to Eq. (23), for (n^) 



as a function of quench time. The transition from KZ scaling to exponential suppression 
is very clear. 

As we saw above, we have G(0) oc / dk k^^^P{k)] it and other moments of P{k) 
play an important role in winding number evolution [TH] . In Fig |3] we take the first 
moment 

tdkP(k) 
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for a given tq. The cutoff k = 2tt corresponds to a length cutoff at 0(^o) in our 
dimensionless units. In Fig. [3] (a) we plot kP{k) from the numerical simulation for 
t > in multiples of ^/tq- As we expect, the power gets pushed into longer and 
longer wavelengths as the field becomes ordered. The scale over which that is ordered 
becomes clearer in the moment I{t) of (b), in which we compare the simulation to 
the analytic behaviour of (15). The expression (15) corresponded to an unbounded 
linear quench without the brakes imposed for \t\ > tq of (25). The comparison to 
simulations is only valid for t > —tq, where it is seen to be good. As it stands the 
Gaussian approximation drives the power to longer wavelengths than the simulation 
but, as we know from elsewhere [TB], the scaling behaviour of the dominant wavelength 
(corresponding to ^~^) still satisfies the same scaling behaviour as for the simulation for 
large annuli. 

We also find that, provided the time is not too late, the diagonalisation of (12) in 
field components is approximately satisfied, but with large noise; we do not present any 
data to that effect. 



3. Simulations of Annular superconductors in three dimensions 

Even thin annuli are not one-dimensional and we should take the annulus width into 
account. Moreover, to compound the problem, even a superconductor that is effectively 
a 2D film cannot be treated entirely as two-dimensional. The three-dimensional nature 
of the flux lines cannot be ignored. In fact, the KZ analysis is known to be incomplete, in 
that the possibility exists of long wavelength modes of the electromagnetic field freezing 
in as the transition is implemented. We have seen elsewhere |8] that this secondary 
mechanism, the Hindmarsh-Rajantie (HR) mechanism [121 [20], does not seem to be 
important for our small systems and we shall not pursue it further. 



3.1. The 2D Gaussian approximation 

For the moment we restrict ourselves to the complex field 0a (x), extended now to the 
Xi —X2 plane, assumed Gaussian, as before. As it stands it is too difficult to get analytic 
results for arbitrary annuli SO, cLS cL first step, we consider a 2D superconducting film of 
infinite extent and we look for the behaviour of the winding number along a circular 
closed loop in the film of circumference C, that encloses a surface S. 

If the topological density [211 ES] in the plane is P3(x), for given field configurations 
0a (x) the phase change 6c along the path can be expressed as the surface integral 

ec = 2Tx I d^x p^ipi) = -27T / d^x p3(x) (29) 

from charge conservation. It then follows from the results of Refs. [211 122] that, for 
Gaussian fields, (n^) satisfies 

{n') = - [ d'x [ d'y (P3(x)p3(y)), (30) 
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Figure 3. (a) We plot kP{k) against k for a variety of i > 0, in multiples of the KZ 
time i> = y/rg for tq = 32 on a lattice of size C = 256. The characteristic wavelength 
increases with increasing time as the field orders itself. Details of the quench are given 
below, (b) We plot the first moment of P{k) from the simulation and compare it 



to its value from (15). The comparison is only valid for t > —tq, where it is good. 
We note that an identifiable scaling length is only established for t > 3i>. This is 
reasonable from the KZ viewpoint since is the minimum time at which the system 
could unfreeze. 
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where x and y are in the plane of S. It is not difficuh to express the density correlation 
in terms of the field correlation /(r) (r = |x|), 

(P3(x)P3(o)) = ^^r-q^). (31) 



Att'^t dr\l — f^ir] 



It follows from (17) that the density correlation has a range 0{^r)- With x outside 5, 
and y inside S, all the contribution to (n^) comes from the vicinity of the boundary 
of 5", rather than the whole area. For large rings this means that, if we removed all 
material except for a strip of width O(^) from the neighbourhood of the contour C we 
would still have the same result. As before, this gives the perimeter rule 

{n')^a'^. (32) 

That is, we have the anticipated random walk in phase along the contour, from which 
we recover the canonical scaling ^ on identifying ^r- 

We can say more. The analysis that led to the expressions (14) and (15) for the 
power spectrum of the fluctuations is dimension independent and only depends on the 
Gaussian nature of the fields. Dimension only arises in the construction (13) for the 
correlation functions, and we shall not repeat the simulations in this case. 

For small annuli the situation is less clear. While the behaviour of (30) suggests 
that the winding number along a small closed path in a 2D superconducting film is 
proportional to the area of the ring (and hence doubles the KZ exponent), it is not clear 
that this behaviour would be preserved if the interior of the closed path were removed 
to make a broad annulus. A test of this is to preserve the perimeter of the annulus 
but to take different aspect ratios (i.e. different areas). A perimeter law would leave 
probabilities unchanged, but an area law would lead to differences. The relevance of this 
is that in an earlier experiment [23] , results suggested a doubling of the KZ exponents - 
albeit with large errors. We shall return to the question of exponent doubling for small 
systems in our conclusions. 



3.2. Simulations 

To test the approximations outlined in the previous sections we have simulated the 
system illustrated in Figure |4| in which the superconducting ring (where the complex 
scalar field is defined) is surrounded on all sides by a box in which the standard U(l) 
vacuum electrodynamics is simulated. Although the superconductor is thereby taken to 
be a thin planar film, the three-dimensional simulation allows for nontrivial correlations 
of the magnetic field. Details are given in [8], to which we refer the reader, and we only 
summarise them here. 

Extending the approach used for the ID ring, we use the gauge-invariant Langevin 
formulation for gauge and scalar fields denoting the lattice order parameter field 
again as and the gauge fields as Ai. 
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Figure 4. The 3D simulation consists of a thin superconducting ring (shown in grey) 
suspended in a box of size 30 x 30 x 5. (we take Sx = 0.5). The ring is of width 5, the 
borders also of thickness 5, while the square central void has size Ax — Ay ~ 10. An 
external magnetic flux density i?oxt is applied, as shown. Periodic boundary conditions 
are used. 



The Hamiltonian for the system is 



H 



jk 



+ 2 ^ Re (j)'^{x)Ui{x)(f){x + 



+ {sxrYl 



e(t) 



{6x) 



2 



(33) 



where Ui{x) = e^^-^^i^^) ^ and the A"*" denote finite differences. Note that is only 
defined in a 2-dimensional superconducting film region; 11 and Ai are coupled in a 
gauge-invariant manner to Gaussian white noise with strength rj. We will take e = 0.1 
and h = 0.1, meaning we are simulating a type II superconductor with k = a/A/c ~ 3.2. 
Periodic boundary conditions are used. The quench protocol is the same as for the ID 



simulations, meaning e{t) is of the form given in Eq. (25). 
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Figure 5. Results of simulations for annuli with inner circumference C in ID and 3D, 
as a function of tq/C^. The ID simulations are represented by open symbols (black), 
the 3D by crosses (red). In the former case we see universal behaviour, and both a 
change in behaviour from KZ scaling to exponential suppression occurs. We see that 
there is complete agreement for large annuli (fast quenches) and different exponential 
damping for small annuli (slow quenches). A possible explanation for the difference 
lies in the way in which field correlations play a role in different dimensions. 



Here, we first sought to see if tlie exponential falloff in trapping probability predicted 
above by ID simulations persisted for the more realistic theory of 2D rings in three 
dimensions. We had in mind the possibility of exponent doubling, suggested earlier 
[23] . In Figure [sj we show the crossover from KZ behaviour to exponential damping as 
a function of the combination rg/C^, for both the ID and 3D systems discussed above. 

Exhibiting universal behaviour of the form (23), the 3D data are represented by 
red crosses superimposed on ID data for a variety of ring circumferences C. If we were 
to plot the same data on log-linear axes we would see that the collapse in probabilities 
is exponential (as we saw in ID in [7]). However, while there is complete agreement for 
large annuli the different exponential damping for small annuli is dimension dependent. 
This is perhaps not surprising because of the differences between the ID expression 
Eq. (16) and the 2D expression (31) on substituting Eq. (20) for /(r), even though the 



annulus is not particularly wide [IH]. As we would have anticipated, this difference is 
insensitive to self-coupling strength, as we have checked (but not displayed). 

Further, in Ref. [H], differing aspect ratios were tried for the central hole for the 
same perimeter, ranging from a square to an elongated 5:1 rectangle, but no deviation 
from perimeter law behaviour ([s]) was detected, leading us to ignore the HR mechanism, 
whose contribution does not obey a perimeter law. 
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4. External fields 

All the analysis above assumed no external fields. Experimentally, this is unrealistic. 
There are nearly always stray fields in the equipment or the environment. Moreover, 
driven symmetry breaking through the application of external bias is an interesting 
question in its own right. The KZ scenario cannot address this, but we shall find that 
we can obtain results in the Gaussian approximation that can be analysed in simulations. 

Experimentally we compensate for unknown stray fields by applying an external 
bias field to cancel them. We know when this has happened when the likelihood of 
seeing spontaneous flux is at a minimum. To do this successfully we need to know the 
likelihood of fiuxoid production in the presence of external fields. In practice it is difficult 
to extend the previous analysis with periodic boundary conditions to include external 
fields and we adopt a simpler straightforward variant of the Gaussian approximation. 



Fiuxoid production in an external field: ID approximation 



Let us now apply a magnetic flux $j to the loop; this breaks the 9 — ?■ —9 symmetry. 
The effect of $/ is to produce a non-zero average winding number (n) = = $//$o 

after the transition, where $o = hc/2e. 

Because of the limitations of one dimension, the electromagnetic field can only be 
introduced classically. The natural extension [15] of our Gaussian probability model 
Gat (6) in Eq. ([T]) for a superconducting loop linked to a magnetic flux $/ is that the 
phase distribution will still be normal with the same variance a'^{N), but with non-zero 
average 6(<l>/) = 27rn($j) = 27r$//$o = 27r0/: 



exp 



(0 - 27r0;)^ 



(34) 



(AT) ' 2a\N) 

With such a distribution, the probability Pm of ending up with a given winding number 
m is given by 

■■7r+2m7r 



Pr. 



d9 G{9](pi). 



7r+2m,7r 



That is 



P±m{N, (pf) 



erf 



(±2m - 20/ + l)7r 



erf 



(±2m 



which, in the presence of a possible stray or residual flux 

1 
2 



>r, we reparametrise as 



Pm{(j){ 



erf 



{ — (j)r — m + 0.5 — (pr — m — 0.5 
erf 



(35) 



(36) 



(37) 



Sd Sd 

where the dependence on quench time and geometry is parameterised by Sd = 
^j2a'^{N) /2t: . Now, j9+m(0/) ^ p_m(0/). Essentially, what really matters is the 
difference m — 0/; for example, Pm(l) = Pm-i(O). This calculation unfortunately does 
not allow us to determine the dependence of s on tq. Nevertheless, we can fit data 
for different apphed magnetic fields at fixed tq to Eq. (37) with only Sd and 0r as free 
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parameters to yield /q. We will test the validity of this ansatz (with 0r = 0) in our 
simulations. 

We have seen that the Gaussian probability is a good approximation for large rings 
and, as before, we extrapolate to small rings, using (36) as it stands. 



4-2. Simulations with an external field 

Where we wish to impose an external field in our simulations, we do so by twisting a 
single plaquette - located at Xq in the x — y (1 — 2) plane. The effect of this, discussed 
in detail in Ref. [8], is to add the term 

i/ext = ciYl (^i(^) + Mx + ai) -Ai{x + a2) - A2{x)) B,^tS^^\x-xo){3i 

X 

to the Hamiltonian, where i?ext is the applied field, which can take arbitrary values 
because we have used the noncompact formulation of the gauge field A. Of course, the 
equations of motion must be modified to take account of this external field term. Results 



for nonzero external field are shown in Figure [6| with fits of the form in Equation (37) 
also shown. Despite the simple assumptions going into this ansatz, it agrees with the 
results of the simulations very well. 



4-3. Experiments with an external field 

We have performed experiments using the setup of Fig. [T]to detect fluxoids produced in 
JVb annuli after a temperature quench in the presence of an external field. Details are 
given in the appendix. It is sufficient for the conclusion of this paper to say that the 
fluxoid probabilities match the simulations well which, in turn, means that they match 
the Gaussian form given in Eq. (37) extremely well. In Fig. [7]we give two examples of 
frequency dependency on the external field, for parameter values very similar to those 
shown in Fig. |6j Unfortunately, we do not have enough experimental data to show how 
pi varies with tq, except that it is in the exponentially suppressed regime. 



5. Conclusions 



In this paper we have explored the ways in which spontaneous fluxoid production in 
annular superconductors can be used to cast light on the Kibble-Zurek (KZ) mechanism 
in systems that are far from equilibrium. In particular, does it represent the dynamics 
well? 

With experiments both past and present in mind, the simple picture is complicated 
by two effects; small size and possible external fields. The KZ scenario is unable to 
address either of these. However, we know that the basic scaling results ^ also provide 
a good test of the Gaussian behaviour of the order parameter field and we have extended 
the Gaussian approximation to those areas which the KZ scenario cannot reach. 

As it stands the Gaussian approximation is only simply analytic for an idealised 
ID annulus, in which we can include an external magnetic field classically. In fact. 
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Figure 6. Simulated probability pi of observing one defect as a function of normalised 
external magnetic flux 0/ for (a) Tq — 180 yielding s — 0.22 ± 0.01 and (b) Tq — 120 
corresponding to s = 0.39 ± 0.01. The curves are a one-parameter fit to s of the form 



given in Eq. ( 37 ) 



we have adopted two variants of Gaussianity which differ very shghtly: one is more 
convenient for addressing periodicity, while the other for including external fields. We 
have performed simulations based on TDGL theory, both for the same idealised ID 
annulus and for a more realistic annulus with finite width embedded in 3D space, in 
which the electromagnetic field is fully dynamical. Only in this latter case can we 
properly include the effects of external fields, necessary so that stray fields can be 
eliminated in experiment but interesting in its own right. 

Our results are surprising. Where comparable, in major regards we find no 
important difference between the ID and 3D simulations. In particular, the likelihood 
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Figure 7. Experimental probability pi of observing onedefect as a function of 



normalised external magnetic flux t/f/.The curves are a fit of the form given in Eq. (37 1. 
Thedashed lines show the best two-parameter fit to Eq.(37): (a) 0r — 0-752 and 
(b) — —0.028 and Sd — 0.378. The integer m is undetermined. These 



Sd = 0.205; 

values are very close to those of Fig. |6] above. Theexperimental data, the simulated 
data and the Gaussian approximationfor the probability are indistinguishable. 



of fluxoid production depends only on the inner circumference of the annulus, and not 
on its shape, all other things being equal. For large rings we reproduce the canonical 
KZ scaling of ([s]). However, the likelihood of seeing fluxoids for small ring size or slow 
quenches is exponentially damped in a dimensionally dependent way. In the presence of 
the explicit symmetry breaking of external fields we find strong agreement between the 
3D simulations and the ID results of the Gaussian approximation. More to the point, 
we show experimental results for JVb annuli that agree totally with both. 

The implications for the KZ scenario are strong. Rather than think in terms of 
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causal horizons as the constraints on growth of correlations, we see field ordering in terms 
of the fastest possible growth of unstable long wavelength modes, as seen both in the 
Gaussian approximation and in simulations. The validity of the Gaussian approximation 
as used here, and hence the scaling behaviour, is predicated upon damped TDGL 
equations describing the non-equilibrium dynamics of these modes, even though the KZ 
picture does not use this language in its adiahatic — )■ impulse — )■ adiahatic progression. 
In another paper in this volume one of us (RR) has shown that a cold tunable condensate 
can pass through the same sequence of freezing - unfreezing, whose causal bounds 
constrain correlation lengths. Although the KZ picture would suggest the formation 
of defects (vortices) in this case the absence of a supporting Gaussian approximation 
mitigates against this. Numerical simulations of this system are underway. 

We conclude with an important issue that really requires a greater analysis. This is 
the empirical observation of scaling behaviour in which the observed exponents initially 
take double the canonical KZ values, as seemed to be the case in our earlier experiment on 
superconducting annuli and as happened in our experiments on annular Josephson 
tunnel junctions |25l [13]. Further, a recent experiment on defect production in linear 
ionic crystals again shows twice the expected exponent [20] ■ We believe there are 
different reasons for all of these. In [23] we provided a possible Gaussian explanation 
for superconductors in terms of small size effects, which explanation was picked up by 
the authors of [26]. Preliminary analyses here suggest something more complicated. 
As we have seen above, the conditions we imposed in [23j for doubling in ID systems 
are not satisfied for superconducting loops. Unless there is an effect due to annulus 
width (which Gaussianity suggests for wide annuli but which we have been unable to 
reproduce in simulations) our 'doubling' of [23], with its large errors in comparison to 
those of [251 IE] , we would now understand as the onset of the shoulder in Fig. |5] as the 
exponential damping takes hold. 

There is the same lack of doubling for short linear superconductors with open 
boundary conditions, where we might have thought the simple explanation had a better 
chance of working. However, the situation is different for short linear systems with open 
boundary conditions for Z2 (rather than U{1)) symmetry breaking, where we seem to 
see doubling of the exponent in our simulations (not presented here). This is more 
appropriate to linear ionic crystals, giving a possible way to understand the empirical 
exponent doubling seen there [26] . 

The doubling of KZ scaling exponents for fluxon production in annular Josephson 
tunnel Junctions is the most complicated of all. On the one hand, it can be understood 
as due to the fabrication of the junction, a proximity effect arising from the (a priori 
undetermined) deposition of non-superconducting metal on the oxide interface. There is 
possibly another cause, due to the presence of a pinning centre for vortices. If we consider 
a short Josephson tunnel junction as an idealisation of a narrow annular junction with 
a strong pinning centre then the doubling effect seen in Z2 symmetry breaking may 
take place here, for the Z symmetry breaking appropriate to the sine-Gordon equation 
in the Josephson angle 6*, whose winding number counts fluxons [27]. This is under 
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consideration. 

Note added: A similar exponential collapse of defect probability to that shown in Fig.|5] 
has been observed in Josephson vortices in annular condensates by Su S-W, Gou S-C, 



Bradley A, Fialko O and Brand J, arXiv: 1302.3304. We thank J Brand for bringing this 
to our attention. 
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Appendix 

For a thin-film loop, which typically has a width, w, much larger than its thickness, the 
intensity of the radial magnetic field. Hp, at the film surface, as derived from Ampere's 
law, is: 

In deriving the above expression we also assumed w smaller than the ring radius 



(or curvature radius in the case of a non-annular loop). Eq. (A.l) indicates that a 
superconducting loop acts as a flux-to-field transformer; if a magnetic sensor is placed 
above (or below) part of the loop, it will thus detect the local radial field. Hp, and so 
the magnetic flux, $e; linked to the loop, and/or the winding number, n. Pannetier 
et a/. [28] successfully used high-sensitivity giant magnetoresistive sensors. However, in 
this context, the most natural magnetic sensors are planar Josephson tunnel junctions. 
More specifically, the critical current of Long Josephson Tunnel Junction (LJTJ) - for 
which the loop itself constitutes one of the superconducting electrodes - is able to resolve 
flux changes well below the flux quantum [22] and promises to be a competitor of the 



magnetic sensors based on quantum interference [30]. Fig. Al shows a superconducting 
ring acting as the base electrode for two L JTJs (the need for two L JTJs will be discussed 
below). Let us stress that the critical current of a LJTJ is sensitive to the surface field. 
Hp, whereas the superconducting loop is only sensitive to the perpendicular field, H±. 



Figs. A2 (a) and (b) show the quantum levels observed when the system is cooled 
through a controlled Normal-Superconducting (NS) transition in the presence of a 
transverse magnetic field, if^, which is incremented by steps corresponding to ~ 0.02 $o 
jumps in the freezing flux, $/ = fioH±Aioop, where Aioop is the loop area. Once H± is 
removed, each trapped flux quantum results in a small change in Hp which, in turn, 
produces a detectable change in the junction zero-field critical current, Ic- The two 
panels of Fig. X2 refer to the same sample quenched in the same field range, but with 



LJTJ 




Figure Al. The magnetic field associated with the trapped fluxoids can be detected 
by one or more Long Josephson Tunnel Junctions, where the base electrode (dark grey) 
is formed by the ring itself. The top electrodes of the two planar Josephson tunnel 
junctions arc in light gray. 
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Figure A2. Winding number, n, vs. the magnetic transverse field, H±, with which the 
loop was field-cooled through its superconducting transition temperature, Tc ~ 9.1 K. 
(a) slow quench: the temperature rate change, dT/dt, during the cooling through the 
critical temperature was — 0.5K/ms; (b) fast quench: dT/dt ^ — 50K/ms. Arbitrary 
horizontal and vertical offsets. It should be noted that H± is switched off during the 
readout. 



different cooling rates, dT/dt: in the first panel it is about — 500K/s, while it is hundred 
times larger for the second panel where a marked overlap of the quantum states is visible. 
This overlap is due to the non-adiabaticity of the thermal quenches and, as discussed 
in the main text, the quantum transition from one level to the next follows a Gaussian 
probability law. In deriving Figs. A2 a) and (b), we assumed that the (unknown) 
residual magnetic field of our setup was small enough to be trapped in the loop as less 
than one fluxoid. Nevertheless, the presence of any larger (static) magnetic field would 
simply result in a horizontal shift of the field axis or, correspondingly, a vertical shift of 
the winding number, n. In other words, we are able to measure the changes in the loop 
quantum state, but not its absolute value. Luckily, as it was argued in the main text, 
the transition from the n-th to the n + 1-th state is virtually independent on n. 



Appendix A. 1. Experimental setup 

Our setup consisted of a cryoprobe inserted vertically in a commercial LHe dewar. The 
cryoprobe was magnetically shielded by means of two concentric Pb cans and a vacuum 
tight cryoperm can surrounding them and immersed in the LHe bath (T ~ 4.2 K). In 
addition, the measurements were carried out in an RF-shielded environment. In the 
experiments, we used high quality all-Niobium LJTJs fabricated on 4.2 x 3mm^ silicon 
substrates using the trilayer technique in which the Josephson junction is realized as 
a window opened in a Si02 insulator layer. More details on the samples electrical 
and geometrical parameters can be found in Ref.[29]. The chip was mounted on one 
side of a massive Cu block. Inside the outer can a He exchange gas with a pressure 
of about 20mbar provided the thermal link between the Cu block and the LHe bath. 
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Figure A3. Quenching time, tq, vs the duration, At, of the laser pulse for three 
samples. The dashed line shows the best allometric fit. 



The chip was heated above the loop critical temperature, Tc ~ 9.1 K, by a laser beam 
pulse transmitted along an optical fiber down to the back side of the chip. The single 
crystalline Si chip absorbed a large fraction of the incident green light, and its very high 
thermal conductance minimized thermal gradients along the chip. After the laser pulse, 
the excess heat diffuses away from the chip through the thermal contact with the Cu 
block and the He gas inside the can; the chip temperature then relaxes down to the 
bath temperature. 

During the transition from the normal to the superconducting states a calibrated 
magnetic field, H±, was applied in the direction perpendicular to the loop plane by 
means of a superconducting cylindrical coil aligned with the loop axis, while the LJTJ 
was electrically isolated: in fact, both the junction voltage and current leads were 
open during the whole thermal cycle. At the end of each cycle, the transverse field 
was removed and - as previously explained - the possible fluxoids are counted by a 
measurement of the LJTJ's (positive and negative) zero-field critical currents. Due 
to the annihilation of a fluxoid-antifluxoid pair, this method works well as long as 
the chances to have two trapped fluxoids are small. Hundreds of thermal cycles were 
carried out for each value of the trapping field, H± ; the field value was increased in steps 
corresponding to about 0.1 — 0.2$o changes in the freezing magnetic flux, $/, linked to 
the loop, until a total flux variation, A$/-, well above one flux quantum was achieved. 
In order to run batches of several thousands of equal thermal cycles with given constant 
parameters, the automation of thermal cycles was implemented by means of a switching 
unit controlled by a GPIB interface; this also allowed for much more robust statistics to 
be achieved. At the end of each thermal quench, the zero-field junction current-voltage 
characteristic is automatically digitally acquired and stored. Later, the output of each 
thermal cycle is analysed algorithmically to yield statistics for the trapped flux quanta. 
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The temperature dependence of the LJTJ's gap voltage was used as an on-chip 
fast thermometer, to monitor the chip temperature during the thermal cycle and 
the quenching time, tq, was inferred from a a well known fitting and extrapolation 



procedure [3T]. Fig. A3 shows that the quenching time, tq, is proportional to the laser 
pulse duration. At, with the proportionality constant almost independent of the sample 
(and of its mounting procedure). Indeed, this method for reading the loop winding 
number by means of a LJTJ, has been already suggested and adopted in a previous 
work[32], but it became reliable and efficient only when the 5-biased[33] junction was 



replaced by an in-line type one[3ll |35l [29], as depicted in Fig. Al 



Appendix A. 2. Fluxoid readout 

The discrete variation, AJc, of the zero-field critical current associated with each fiux 
quantum trapped in the loop is[ 



/\h = gp^, (A.2) 
where is a dimensionless parameter of the order of unity that depends on some 



geometrical details [29]. With a loop inductance Lioop ~ 100 pH, Eq.(A.2) provides Ic 
jumps of the order of several microamperes in a fiux range of hundreds of fiux quanta|30] . 
The discreteness of the critical current values is wider, when the loop is realized with 
a top junction electrode that, with respect to the bottom electrode, typically, has a 
smaller width and a larger thickness, which both result in a larger gain factor, g 



The validity of Eq.(A.2) has been supported by measurements on a number of devices 
based on narrow loops with annular and rectangular geometries, namely, rings with mean 
radius i? = 53 /im, width w = Q /im and inductance Lioop ~ 240 pH and rectangular loops 
with mean base h = 150 fim, mean height h = 30 /im, width w = 6 fim and inductance 
Lioop ~ 90 pH. 

Appendix A. 3. Common mode 



In the common mode configuration, depicted in Fig. XT, the two LJTJs share the same 
doubly connected base electrode and are series biased. Any circulating current, Idr, 
modulates oppositely the zero-field critical currents. Id and Ic2, of the two LJTJs; 
then, their difference, 6Ic = Ic2 ~ hi, changes twice faster with the winding number, 
n. Further, the rms current noise on 5Ic is only a/2 times larger than that on a single 
critical current, meaning that in the common mode configuration the signal-to- noise 
ratio is enhanced by a factor v^- However, the main advantage of such configuration is 
that by reading the two critical currents simultaneously and looking at the correlation 
between their changes, it is possible to exclude unwanted events from the statistical 
analysis. These are related to the trapping of Abrikosov vortices found at the pinning 
centers of the superconducting electrodes (such as normal inclusions, dislocations and 
grain boundaries). 
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Figure A4. (a) Zero-field critical currents Id (dark points) and /c2 (light points) 
as a function of the cycle number i. (b) The same for Sic = Ic2 ~ Ici- (c) Relative 
variation of the freezing flux, $j. All plots have an arbitrary vertical offset. 



When using two series-biased LJTJs, care should be taken to ensure that their 
critical currents are sufficiently different, otherwise the switching to non-zero voltages 
of the weaker junction might prematurely drive the other to the finite voltage state. 
One further reduction by a \/2 factor of the signal-to-noise ratio is achieved by taking 
the combination, -^(^2 ~ -^c2 ~ -^ci + -^ci' where 2 the absolute values of, 

respectively, the positive and negative critical currents of the two LJTJs. 

Appendix A. 4- Measuring the trapping frequency 

Fig. |A4K a) shows the zero-field critical currents. Id and Ic2, of a two-junction annular 
sample during a statistical batch in which each value of the cooling field, H±, was 
maintained for about thousand thermal cycles. Each cycle is associated with a 
progressive integer i, so that the index i grows linearly with time. Since each cycle lasts 
slightly less that 4 seconds, the total batch, consisting of 15000 cycles, required about 
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twelve hours. The reduction of the critical currents observed during the measurements 
corresponds to a small decrease (drift) of the bath temperature. The middle panel of 
Fig. A4 plots the difference 6Ic = Ic2 — hi, we observe that the current jumps have 
doubled, the signal-to-noise has been somehow lowered and, at the same time, the effect 
of the temperature drift has been counterbalanced. As Fig. |A4](c) indicated, the cycle 
number i shows the step-like dependence of the freezing flux, as time goes by. Let us 
call the residual magnetic flux; then, to reproduce the zero-fiux condition during the 
quench, a freezing flux $/ = — $r needs to be applied. Since is, in general, unknown, 
it is necessary to span $ j over a range at least as large as one flux quantum. 

From the statistical analysis of data this obtained - such as that shown in Fig.|A4](b) 
- it is possible to determine the frequency distributions, fm- In order to compare loops 
with different effective area, the trapping frequency, fm, is plotted against the normalized 
freezing magnetic flux, 0/ = ^f/^o, rather than the freezing field, H±. Two examples 
of the computed frequency distributions were given in the main text in Figs. [7|^a) and 
(b); the statistical data are very nicely fit by a Gaussian behavior, as in Eq.(|37]). 
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